The theory of optical dispersive shock waves in photorefractive media 
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The theory of optical dispersive shocks generated in propagation of light beams through photore- 
fractive media is developed. Full one-dimensional analytical theory based on the Whitham modu- 
lation approach is given for the simplest case of sharp step-like initial discontinuity in a beam with 
one-dimensional strip-like geometry. This approach is confirmed by numerical simulations which are 
extended also to beams with cylindrical symmetry. The theory explains recent experiments where 
such dispersive shock waves have been observed. 
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I. INTRODUCTION 



Study of optical solitons is a large area of modern research which is important both scientifically and for potential 
applications (see, e.g., P, 0). Different kinds of solitons have already been observed in various nonlinear optical 
media and their behavior has been explained in the frameworks of such mathematical models as nonlinear Schrodinger 
(NLS) and generalized nonlinear Schrodinger (GNLS) equations for different dimensions and geometries, so that one 
can consider the properties of single solitons as well enough understood. 

However, there are situations when many solitons are generated so that they comprise a dense soliton lattice. In such 
situations, it is impossible to neglect interactions between solitons and one has to consider evolution of the structure as 
a whole rather than to trace evolution of each soliton separately. Usually, such soliton structures appear as a result of 
wave breaking of a large enough initial pulse or large disturbance along constant background. Hence, such structures 
can be considered as a dispersive counterparts of shock waves well known in physics of compressible viscous fluids 
(see, e.g., Q). In a viscous fluid, the shock can be represented as a narrow region within which strong dissipation 
processes take place. In optics, on the contrary, dissipation effects can be neglected compared with dispersion ones 
and the shock discontinuity resolves into an expanding region filled with nonlinear oscillations. Such dispersive shock 
waves are known as tidal bores in rivers [3] and have been also observed in some other physical systems including 
collisionless plasma [5j and Bose-Einstein condensate [6| . Experiments on such dispersive shock waves production in 
optics have been recently reported in @, HI- Motivated by these experiments, we shall consider here the theory of 
dispersive shock waves in photorefractive media. 

Since the number of interacting solitons in dispersive shocks is usually much greater than unity and these solitons 
are spatially ranked in amplitude, such a dispersive shock can be represented as a modulated periodic wave with 
parameters changing a little in one transverse or longitudinal period of envelope amplitude of the electromagnetic 
wave. A slow change of the parameters of the envelope amplitude is governed to leading order by the Whitham 
modulation equations obtained by averaging conservation laws over the family of nonlinear periodic solutions or by 
the application of the averaged variational principle (see, e.g., [3, 13, [lOj ) . For the one-dimensional NLS equation, the 
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Whitham equations were derived in [ll|, 1 1211 (see als o 11011) and the mathematical theory of dispersive shock waves 
for the defocusing case was developed in jl3l . Q, [ll| Il6l. \VA llR [l9|. It was applied to the propagation of signals in 
optical fibers in [20| and in Bose- Einstein condensates in [y, [2l|. It should be mentioned that for the case of the ID 
NLS equation, the presence of an integrable structure has important consequences for the modulation (Whitham) 
system, namely, the latter can be represented in a diagonal (Ricmann) form, which dramatically simplifies further 
analysis. The method of obtaining the Whitham equations in this form is based on the Inverse Scattering Transform 
(1ST) applied to the NLS equation (HI, [12]. However, in the case of the GNLS equation, the 1ST method cannot 
be used anymore, and the diagonal structure of the Whitham system is not available. Nevertheless, it was shown in 
[2^-[24| that in this case, the main characteristics of the dispersive shock wave still can be found by using some general 
properties of the Whitham equations which remain present even in non-integrable case. Here we shall use this latter 
method for derivation of parameters of one-dimensional dispersive shock waves generated in photorefractive crystals 
and shall confirm our analytical results by numerical simulations, which also provide a more detailed information in 
the cases when the analytical approach is not yet developed (say in 2D). 



II. MAIN EQUATIONS 

Photorefractive optical solitons were first observed in the experiment (25[ and in the experiments @, [1] the formation 
of dispersive shock waves has been observed in spatial evolution of light beams propagating through self-defocusing 
photorefractive crystals, so that beam non-uniformities give rise to breaking singularities and their resolution through 
dispersive shocks. As is known, propagation of such stationary beams is described by the equation 

^ + -La^ + -MM 2 )^ = o, (i) 

oz 2k n v ' 

where tp is envelope field strength of electromagnetic wave with wave number ko = 2irno/\, z is the coordinate along 
the beam, x,y are transverse coordinates, r = (x,y), A± = d 2 /d 2 x + d 2 /d 2 y is transverse Laplacian, n is a linear 
refractive index, and in photo-refractive medium we have 



5n = -\nlr^E p -^—, (2) 

2 P + Pd 



where E p is applied electric field, r 33 electro-optical index, p = IV'I 2 , and pd is a saturation parameter. 
For mathematical convenience, wc introduce non-dimensional variables 



z = ^knlr 33 E p ^j^j z, x = kn ^j ^r 33 E p {j^j x i V = kn o ]j 7} r 33 E P (j^j V: ^ = Va^, ( 3 ) 

where p c is a characteristic value of optical intensity (its concrete definition depends on the problem under consider- 
ation; for instance, it can the background intensity), so that Eq. takes the form of GNLS equation 



; aI + 2 A ^-IT^# 



^ + ±A^-^i-^ = 0, (4) 



where 7 = p c /pd and tildes are omitted. If saturation effect is negligibly small (7|V'| 2 <C 1), then this equation reduces 
to the usual NLS equation 

i|£ + 1a.lV - = 0. (5) 
oz 2 

It is convenient to represent these equations in a fluid dynamics type form by means of the substitution 



V>0, z) = yfp cxp (i J u(r, z)drj , 



(6) 



so that they are transformed to 

p z + V_l(/ju) = 0, 



u z + (uV x ) u + V j_f(p) - V x 



Axp (V ±P ) 2 1 „ (7) 



Ap 8p 2 



= 0, 
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where 

ftp) = — — for the GNLS equation Q (8) 
1 + IP 

and 

f(p) = p for the NLS equation © . (9) 

The light intensity p in the hydrodynamic interpretation has a meaning of a density of a "fluid" and Eqs. I|8|). ([9]) can 
be viewed as "equations of state" for such a fluid. The function u(r, z) is a local value of the wave vector component 
transverse to the direction of the light beam; in hydrodynamic representation it has a meaning of the "flow velocity" . 
The variable z plays the role of time so it is natural to describe the deformations of the light beam in evolutionary 
terms. Obviously, if initial distribution does not depend on one of the transverse coordinates (say y), then transverse 
differential vector operators reduce to usual derivatives (Vj_ = d/dx, Aj_ = d 2 /dx 2 ). 

Evolution, according to ([7]) of an initial distribution, specified at z = 0, typically leads to wave breaking and 
formation of dispersive shock waves. One can distinguish the following typical cases: 

• generation of dispersive shocks in the evolution of a bright strip hump above a uniform (background) intensity 
distribution, 

• generation of sequences of solitons from a strip "hole" in the light intensity, 

• generation of dispersive shocks in the evolution of a bright cylindrically symmetrical hump above a uniform 
intensity distribution. 

In ID geometry such humps can be modeled qualitatively by step-like pulses with sharp boundaries, and these models 
are convenient for analytical considerations. As was shown in [6| for the NLS equation case with 7 = 0, this model 
agrees quite well with numerical simulations of 2D dynamics. Therefore we shall start here with these idealized 
models. 



III. ANALYTICAL THEORY OF ONE-DIMENSIONAL DISPERSIVE SHOCKS GENERATED IN 

DECAY OF A STEP-LIKE INITIAL DISTRIBUTION 



We shall start with analytical treatment of shocks described by ID equation 

1 



np z + ^ xx -fM 2 )iP = o, (10) 



or, in a fluid dynamics form, by the system 



Pz + {pu) x = 0, 
uu x + —p x + -—- — = 0, 



dp r V 8 /° 2 4 P 

where the nonlinear refraction function f(p) is given by Eq. ([5]) or ©. The systems of the type (fTI]) are often referred 
to as dispersive hydrodynamics systems. 

We consider initial distributions of the intensity and transverse wave vector in the form 

p(ar,0) = |{° u(x,0)=0, (12) 

that is we assume that the initial velocity u{x, 0) is equal to zero everywhere which means that the initial beam enters 
the photo-refractive medium at z — without any focusing. For the sake of definiteness we assume also that po > 1. 

At the initial stage of evolution, linear waves are generated which propagate according to the dispersion law obtained 
by means of linearization of Eqs. (JTTJ) about the uniform state p — po, u = u$ (we keep here a nonzero value of uq 
for future convenience), that is p = p$ + p\ exp[i(kx — luz)], u = uq + u\ exp[i(kx — u)z)\, where p\, u% <C 1. Then a 
simple calculation yields 



uj = aj (po,u ,k) = ku ± kJ P ° — Y2+-T- (I 3 ) 

V (i + 1P0) 4 
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Note that lo" (k) > 0, which implies appearance of dark solitons in full nonlinear solutions. But before consideration 
of such solutions, we shall discuss a nonlinear stage of evolution in dispcrsionless approximation when one can neglect 
the higher order terms in the system (fTTjl . While in the case of general smooth initial data this stage of evolution 
is responsible for the formation of breaking singularities in the solution, its consideration also provides important 
insights into the nonlinear dissipationless dispersive dynamics of discontinuous disturbances of the type 1|12[) even 
beyond the breaking point. 



A. Dispersionless approximation 

In dispersionless approximation, the system ([TTj) reduces to standard equations of compressible fluid dynamics 

pz + {f>u)x = 0, 



, „ ( 14 ) 

u z + uu x + } (p)p x = 0. 

Because of the bi-directional nature of this system, generally, an initial step (|12p resolves into a combination of 
two waves propagating in opposite directions. One of these waves represents a rarefaction wave with clear physical 
meaning, but the other one leads to a multi-valued dependence of the intensity p(x, z) and transverse wave number 
(associated flow velocity) u(x, z) on the rr-coordinate. Nevertheless, this formal global solution sheds some light on 
the structure of the actual physical solution and some its elements will be used later, therefore we shall consider it 
here. To this end we cast the system (TT4|) into a diagonal form (see, for instance, p|, [l(|) by introduction of new 
variables — Riemann invariants 

2 

r± = u ± — arctan -i/tp, (15) 

so that it takes the form 

dr+ dr+ , . 

where the characteristic velocities V± are expressed in terms of the hydrodynamic variables p, u by the relationships 

V ± =u±-^-. (17) 
1 +1P 

When 7 — > we have r± — u ± 2^/p, V± — u ± ^fp, i.e. the usual expressions for the dispersionless limit of the 
defocusing NLS equation (the shallow- water system — see, for instance, [l3|). 

Since in the case of the step-like initial conditions the variables r± must depend on a self-similar variable C — x/ z 
alone, the equations Eqs. (fll)]) reduce to (V± — Q{dr±/ dQ = and we arrive at the so-called simple-wave solutions: 

VP x 2 / — _ o 



1+jpz ^7 



arctan-y/Tp = r_ = constant, (18) 



J~p x 2 n 

u = — , u-\ — arctan ^Hp = r, = constant. (19) 

1 + 7P z 

The constants here are chosen from the continuity conditions at the points where the simple waves enter the regions 
of constant intensities. Since the left-propagating rarefaction wave described by (|19p matches with the external flow 
p = po, u = (see Fig. la) we have r° = -j= arctan y/^fpo and, correspondingly 

2 

u = (arctan -v/TPo — arctan y^jp) . (20) 

Now, substituting this into the first equation (fH)]) we get 

i~p 2 x 

— 1 (arctan J^fp — arctan \ffpa) = > (21) 

1 + 7P yfl z 
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which determines implicitly the intensity p as a function of x/z in the rarefaction wave. For x < x^ we have 
p = p Q = constant so x = x^ is the point of weak discontinuity which must propagate with sound velocity (see, for 
instance [26]) which in our case is 



Cs{p) = 



Vp ; 

1 + 7P ! 



(22) 



Indeed, substituting p — po into (f2Tj) we get x^ j z = —c s (po). As a matter of fact, the speeds of propagation of weak 
discontinuities in the photorefractive system agree with the group speeds determined by the long wavelength limit 
k — > in the linear dispersion relation (|13[) . 

Next, for x > x^ we have p = 1, it = (see Fig. la) and this does not agree with the relationship ([20| in the 
constructed left-propagating rarefaction wave solution. Hence, we have to introduce some intermediate distribution 



p(x/z) = p = constant, u{x/ z) 



constant 



(23) 



which matches with the rarefaction wave at some x = x^ . Now, to connect the intermediate distribution 
with p = 1, u = downstream, we have to use the right-propagating simple wave solution (|18[) where the constant 
r° = arctan ^/7- Hence we get 



v/7 



(arctan ^TP — arctan ^y) 



and 



Vp 



(arctan ^ — arctan Vtp) 



x 
z 



Equations (f2"0"|) and 



1 + 7P V"i 

at p = p~ must give it = u~; hence they yield the equation 



arctan V IP 



l . , , v 

— (arctan ^/7po + arctan ^7) 



which determines the parameter p : 



P = 



VI + 7Po - 1 + - 1) 



7VAJ - (V1+7P0 - l)(VT+7 - 1) 



When p is known, the parameter it is found from the equation (|24[) . 

it - = f arctan \/7P~ 

\/7 v 



arctan 1^/7 



(24) 



(25) 



(26) 



(27) 



(28) 



The "internal" end points xf and x 2 are found by substituting the intermediate values p , u into the similarity 
solutions UHl), (fTU), 



' P" 



1 



1P~ 



1 



■7p- 



(29) 



These points correspond to the weak discontinuities which propagate with sound velocities c s (p~~ ) in opposite directions 
in the reference frame associated with the uniform flow u~ . The whole structure of intensity distribution is shown in 
Fig. la. It has the region x^ < x < x\ with the three- valued intensity, corresponding to the formal solution (fTSj) . 
which is obviously non-physical and its appearance serves as an indication that an oscillating dispersive shock wave is 
generated in the region of transition from p = p~ , it = it~ to p + = 1, u + = 0. The arising physical structure is shown 
schematically in Fig. lb. Importantly, the boundaries x\^ of the oscillatory zone by no means coincide with those in 
the formal three- valued dispersionless solution. It is remarkable, however, that in spite of such radical qualitative and 
quantitative change of the flow, the values of p~ and u~ themselves turn out to be still determined by the previous 
equations (j2"?| and (f2"B"l) . This is a consequence of the dispersive shock jump condition which requires that the values 
of the Riemann invariant ?*_ = u — (2/ ^/j) arctan ^7P & t both end points of the dispersive shock wave must be equal 
to each other: 



r -\ x - = r-\ x + , 



(30) 
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p = p 





FIG. 1: Decay of the initial discontinuity of light intensity in a beam propagating through a photorefractive crystal, (a) 
Dispersionless approximation with non-physical region of multi-valued intensity, (b) Schematic picture of formation of dispersive 
shock due to interplay of dispersive and nonlinear effects. The values of x\ and xf are the same for (a) and (b) while the 
are different. 



which gives at once Eq. (|28p . Since the rarefaction wave, even in the presence of dispersion, is still described with good 
accuracy by the dispersionless approximation (see [27| , |28j for the general linear asymptotic analysis of the dispersive 
resolution of the weak discontinuities at the edges of the rarefaction wave), we deduce that Eq. (f27|) obtained in the 
framework of the dispersionless fluid dynamics also remains valid. One should emphasize that, although all obtained 
relationships, strictly speaking, hold only asymptotically for sufficiently large "times" z, as we shall see from the 
direct numerical solution, they hold with good accuracy for rather moderate z. The dispersive jump condition of the 
type pop was proposed for the first time in [28j where it was based on intuitive physical reasoning and the results of 
numerical simulations of collisionlcss plasma flows. A consistent mathematical derivation of this condition along with 
some important restrictions to its applicability was given in the framework of the Whitham theory in [22I |24| . 

As was mentioned, the end points of the oscillatory region of the dispersive shock in Fig. lb do not coincide with 
the end points of the three-valued region in Fig. la. Indeed, this oscillatory zone arises due to interplay of dispersion 
and nonlinear effects and has a structure similar to that observed in the much studied integrable defocusing NLS 
equation case (see [l3| -[ll|). Namely, near the leading edge x^ the wave transforms into a vanishing amplitude 
linear wavepacket and at the trailing edge it converts into a dark soliton. Hence, the end point of the oscillatory 
zone X2 must move with the group velocity of linear waves c g = dojo/dk calculated for some non-zero value of 
k = k + in contrast to the dispersionless approximation corresponding to k — > (in addition to vanishing amplitude 
of oscillations a — > 0). The end point moves with the corresponding soliton velocity which also has nothing to do 
with the dispersionless limit (note, that in the soliton limit k — > but the amplitude a = a~ remains finite). Thus, 
our task is to determine the main quantitative characteristics of the oscillatory region of the dispersive shock — the 
velocities of its end points as well as the amplitude a~ of the trailing soliton at x — x^ and the wave number fc + at 
the leading edge point x = x\ . 

One can observe that the oscillatory structure of the dispersive shock wave is characterized by two different spatial 
scales: the intensity oscillates very fast inside the shock but the parameters of the fast oscillations change little in one 
wavelength in x-direction and in one period along the beam z-axis. This suggests that the oscillatory dispersive shock 
can be represented as a slowly modulated nonlinear periodic wave and, hence, we can apply the Whitham modulation 
theory [3] to its description. In the Whitham approach, the original equation containing higher order x-derivatives 
is averaged over the family of nonlinear periodic traveling wave solutions. As a result, one obtains a system of the 
first order nonlinear partial differential equations of hydrodynamic type (i.e. linear with respect to first derivatives) 
governing the slow evolution of modulations. The modulation system does not contain any parameters of the length 
dimension, so it allows one to introduce the edges x^(z) of the dispersive shock wave in a mathematically consistent 
way, as characteristics where matching of the "internal" (modulation) and "external" (dispersionless fluid dynamics) 
solutions occurs. Of course, strictly speaking the averaged description is valid only when the ratio of the typical 
wavelength to the width of the oscillatory zone is small. For our case of the decay of an initial discontinuity this 
corresponds to a "long-time" asymptotic behaviour, z> 1. However, as we shall see from the comparison with direct 
numerical solution, the results of the modulation approach turn out to be valid even for rather moderate values of z. 

The modulation approach to the description of dispersive shock waves was realized for the first time by Gurevich 
and Pitaevskii [27] in the framework of the Korteweg - de Vries equation. To put this approach into practice for the 
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light beam deformations in a photorefractive medium, we first have to study periodic solutions of the equations (111)) . 



B. Periodic waves and solitons in photorefractive crystals 

The traveling wave solution of the system (|11[) is obtained by the substitution p — p(0), u = u(9), where 9 = x — cz 
is the phase and c = constant is the phase velocity. As a result, we obtain by integrating the first equation (jlip 

u = c+-, 31) 
P 

where A is an arbitrary constant. Substituting (|3ip into the second equation (|lip and performing one integration 
with respect to 6 we obtain an ordinary differential equation of the second order, 

1 ( dp\ 2 ld 2 p 2 2 A 2 

P~ P f(P)-Bp ~ — , (32) 



J9 J AdO 2 ' 1 2 ' 

where B is another constant of integration. We shall seek its integral in the form 



dp \ ^ f 

jq) =a 1 p f(p)dp + a 2 p 2 + a 3 p + a i , (33) 

where a%, a 2 , 03, and 04 are the constant coefficients to be found. Substituting (|33[) into (|32[) we find, with the account 
of the specific dependence f(p), the eventual form of the sought integral, 

^) =-^ln(l + 7/9)+ (a 2 + ^ p 2 + a 3 p + ai = Q(p). (34) 

Here a 2 , a 3 and 04 are arbitrary constants two of which are connected with A and B by the relations 

a 2 = 8B, a 4 = -AA 2 , (35) 

and 03 is an additional constant so that (|34[1 is indeed the first integral of Eq. (|32|) . We denote the roots of the 
equation Q{p) — as e\ < e 2 < e 3 . Then the density oscillations in the traveling wave occur between e\ and e 2 . The 
amplitude of the wave is then given by a = e% — e\. The small-amplitude linear wave configuration corresponds to 
ei — > e 2 while for solitons we have e 2 — e 3 . By imposing the periodicity condition p(9) = p(9 + 27r/fc) we find the 
wave number k of the traveling wave in the form of the integral 



n {L Void) 



(36) 



While the equation (l34|) cannot be integrated in closed form, it is not difficult to find the relationships characterizing 
its special solution in the form of a dark soliton. For this solution we must have the following boundary conditions 
satisfied at infinity: 

p^ p b , u b , dp/dO -► 0, d 2 p/dd 2 -> for \d\ -► 00, (37) 

plus the condition dp/ d6 — at p = p m < pi,, where p m is the value of the "density" in the minimum of the dark 
soliton and pb is the "background" intensity. Applying these conditions to (|3ip . (|34|) we obtain, after simple algebra, 
the expressions for the coefficients in (|34"]) for the soliton configuration, 



a 2 = -— 4(u fc - c) 2 , 

1 + IPb 

s 8p b 4(u b -c) 2 {p 2 m + p 2 ) (38) 

a 3 = — lnl + p b — H , 

7 7(! + 7Pfe) Pb 

a 4 = -4(ub - c) 2 p\. 

The curves Q(p) in a "soliton configuration" for several values of 7 are shown in Fig. 2. The condition that in the 
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soliton limit pb is a double zero of the function Q(p), that is dQ(p)/dp — at p — pb, yields the relationship between 
the soliton velocity c and the amplitude a = pb — p m for given pb, Ub'. 



(c - u b ) 2 



2p„ 



j a 



7a 1 - 



iPb 



1 



IPr, 



1 



IPb 



(39) 



The dependence of the soliton velocity on the saturation parameter 7 is shown in Fig. 3. 

For future analysis it is important to introduce one more parameter — the inverse half-width k of the soliton — using 
the exponential decay of the intensity pb — p as |0| — » 00: 



00. 



Pb — p oc exp(— k|0|), \0\ 
tat 

p' = p b -p and find {dp'/dO) 2 = {l/2){d 2 Q/dp 2 ) Pb {p') 2 = K 2 (p') 2 ; hence 

8p m + ^lPb{pb + p m ) 8p„ 



(40) 



To find the relationship between /c and other parameters we take the series expansion of Q(p) for small values of 



( 1 d 2 Q 


\ 1/2 


y2 dp 2 


Pb) 



In 



1 + IPb 



l(pb- Prn){l + 1 Pm) 2 J 2 (Pb ~ Pm) 2 1+7/°' 

The dependence of n on 7 is shown in Fig. 4. 

The profile of the intensity p(9) is determined by the integral (see Eq. I|34p) 

dp 



1/2 



(41) 



Vow' 



(42) 
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2 4 6 8 10 1 

FIG. 4: The plot of inverse half- width k of photorefractive soliton as a function of saturation parameter 7. The other parameters 
are: p b = 1, p m = 0.2. 




FIG. 5: Profiles of the intensity in photorefractive solitons for values of 7 = 0, 1, 2. The other parameters are: pb = 1, p m = 0.2. 



where it is assumed that the intensity p takes the minimal value p = p m at 9 = which determines the integration 
constant. The wave form of a dark soliton for different values of the parameter 7 is shown in Fig. 5. 
For 7 <C 1 we have asymptotic expansions (for simplicity we take u& = 0) 



C = V^T (l - ~(2p 6 + p m )) + 0(j 2 ). 



K = 2^p b - p m l-^(3p b +Pm) +C>(7 2 ): 



and for 7 ^> 1 other expansions 



c = v/2[/? m (ln(p b /p m ) - 1) + /4/Vb] + -2n 



(43) 
(44) 

(45) 



K = 2 VPb ~ 2p m \n{p b /p m ) - p 2 m /p b + ^ _ 3 

(P6 - Pm)7 

One can see that the leading terms in (|43|) and (|44]) agree with the well-known dependencies for dark solitons of 
the NLS equation 

The particular case of soliton solution with p rn — 0, u b = (hence c = 0) in photorefractive media has been found 
in [H. 



10 



C. Dispersive shock wave 



The general periodic solution of the photorefractive equation depends on the fast phase variable 6 and is character- 
ized by four parameters e\, e 2 , e 3 , c, where e_y, j = 1, 2, 3 are the zeroes of the function Q(p) (|34p . which determine the 
profile of the intensity, and c is the phase velocity. In a modulated wave, these four parameters become slow variables 
of x and z. In the Whitham theory [3!], it is postulated that this slow evolution (modulation) ej(x, z), c(x,z) can be 
found from the conservation laws of the dispersive equation averaged over fast oscillations with respect to the phase 
variable 9. An additional modulation equation naturally arises as the wave number conservation law k z + lo x = 
and essentially represents a condition of the existence of a slowly modulated periodic wave (see, for instance, 0] ) . 
Several averaging procedures have been proposed yielding equivalent results for various physical systems (see [9() so 
the Whitham modulation theory can be now considered as quite well established. As a result, using the original pro- 
cedure of averaging conservation laws, the Whitham system for the GNLS equation can be obtained in the following 
general form 

(Pi(ei,e 2 ,e 3 ,c)) z + {Qi(ei,e 2 ,e 3 ,c)) x = 0, z = 1, 2, 3, (47) 



(fc(ei,e 2 ,e 3 ,c)) 2 + (uj(e 1: e 2 ,e 3 ,c)) x = 0, w = kc. (48) 

Here Pi = p, P 2 = u, P 3 = pu are the conserved "densities" of the GNLS equation {7} and Qi, i = 1,2,3, are the 
corresponding "fluxes" . The averaging is performed over the periodic family (f3Tj) , ([34]) according to 

/(ei,e 2 ,e 3 ,c) = - / == dp. (49) 

n J vQ(p) 

ei 

Now the system (|47p . (|4"51) is, in principle, completely defined. 

The modulation system (|4"T|) , (|4"5|) being the system of hydrodynamic type can be hyperbolic (real characteristic 
velocities - modulationally stable case) or elliptic (complex characteristic velocities - modulationally unstable case). 
It is known very well (see [TTj] , [l^], [Hi) that for the defocusing NLS equation, which is an integrable particular case 
of the the GNLS equation (j4|), the modulation system is strictly hyperbolic. Our numerical simulations show that 
traveling waves in the GNLS equation are modulationally stable and this suggests that the corresponding Whitham 
system is hyperbolic as well. So, in what follows, we shall assume hyperbolicity of the Whitham system, which will 
allow us to use some arguments of classical characteristics theory [3, l26l. l30l] . 

Now, to describe analytically the dispersive shock wave as a whole, we have to solve four modulation equations (|47[) . 
(|48[) for the slowly varying parameters e±, e 2l e 3 and c of the periodic solution. These equations must be equipped with 
special matching conditions to guarantee continuity of the mean flow at the free boundaries x ± (z) defining the edges 
of the dispersive shock wave. In view of the numerically established qualitative spatial structure of the photorefractive 
dispersive shock wave (see Fig. lb) we require that 

at x — x + (z) : a = 0, ~p = p + , u = u + , (50) 



at x = x (z) : k = , p = p , u = u , (51) 

where x + = x J (from now on we shall omit the subscript 2 in x^ and x^). The dependencies of p, u, k, a on 
ei,e 2l e 3 ,c are defined by (|49|) and formulae of Section III.B, and the pairs (p~ , u~) and (p + , u + ) represent the 
solution of the dispersionless approximation (fl4|) evaluated at the trailing and leading edges of the dispersive shock 
wave respectively. The edges x ± (z) of the dispersive shock wave represent free boundaries defined by the kinematic 
boundary conditions with clear physical meaning explained in Section III. A: 

^j— = c g (p + ,u + ,k + ) , < l?— = c so i{p-,u-,a-), (52) 

where c g (p + ,u + ,k) = du^/dk is the group velocity of the linear wave packet with the dominant wavenumber k 
propagating against the hydrodynamic background p + , u + (see (fTB")) for the linear dispersion relation lo = loo{po, u 0jk)) 
and c so i (p~ ,u~ ,a) is the velocity of the dark soliton with amplitude a propagating against the background p~ , u~ 
(see (|3"9")l for the dependence of the soliton velocity on its amplitude). Of course, the values of the wavenumber k + at 
the leading edge and the amplitude a" of the trailing dark soliton are both to be determined, so the determination 
of the edges (z) represents a part of this nonlinear boundary value problem. 
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Following the pioneering work of Gurevich and Pitaevskii [27| on the dispersive shock wave description in the 
framework of the KdV equation, the effective methods for treatment of such problems have been developed for the 
whole class of evolution equations which share with the KdV equation the unique property of complete integrability 
(see, e.g., [HI)- On the level of the Whitham equations, one of the manifestations of integrability is the presence of 
the full system of Riemann invariants, an event generally highly unlikely for the systems of hydrodynamic type with 
number of equations exceeding two. In particular, the NLS equation ([5|) be longs t o this class, and the corresponding 
theory of dispersive shock formation was developed in the papers 



the description of shocks in nonlinear optics [2fJ and Bose-Einstein condensates [g, 

equation (|4]) is not completely integrable and therefore the methods based on the presence of rich underlying algebraic 



181. Il9j and successfully applied to 
2l|. However, the photorefractive 



structure of such equations cannot be applied here. Nevertheless, as was shown in [22[-[2J], the main quantitative 
characteristics of the dispersive shock wave can be derived using the general properties of the Whitham equations (|47[) . 
([4"5| reflecting their origin as certain averages, and here we shall apply this method to the description of dispersive 
shock waves in photorefractive media. To be specific, we shall be interested in the locations of the edges of the 
dispersive shock wave and in the amplitude of the largest (deepest) soliton at the trailing edge, the parameters that 
are usually observed in experiment. 

The method of Refs. [23|-[23|. which will be used below, is formulated most conveniently in terms of the physical 
modulation parameters ~p, u, k 7 a appearing in the matching conditions (|5ip . (|50[) . The key of the method lies in the 
fact that the modulation system (j4"Tj) , (|4"5)) dramatically simplifies in the cases (a = 0, k ^ 0) and (k = 0, a ^ 0) 



corresponding to the limiting wave regimes realized at the boundaries of the dispersive shock wave. 



1. Leading edge 

At the leading edge x = x + (z) the amplitude of oscillations vanishes, a — 0. Since the Whitham averaging procedure 
remains valid for the case a = (averaging over the periodic family with vanishing amplitude), then we conclude that 
the Whitham system must admit an exact reduction at a — and, therefore, the system of four Whitham equations 
must reduce here to only three equations. Now, if the amplitude of oscillations v anishes, then the average of a function 
of the oscillating variable equals to the same function of the averaged variable: F(p, u) = F(~p,u). Thus, when a = 
the Whitham system must agree with the dispersionless approximation (I14|) describing large-scale non-oscillating 
flows, i.e. the modulation equations for p, u, a reduce to 

a = 0, p z + (pu) x = 0, u z + uu x + f'(~p)~p x = 0. (53) 

We note that this reduction of the Whitham equations is also consistent with the matching condition (|50| at the 
leading edge of the dispersive shock wave where a = and which requires that the solution of the Whitham equations 
must match with the solution of the equations of the dispersionless approximation. Of course, equations (|53p can be 
derived directly from the modulation equations (|4T[) by passing in them to the limit a = ei — e\ — > (see, for instance, 
[33l ] for the corresponding calculation in the context of fully nonlinear shallow- water waves), however, validity of (|53p 
appears to be obvious from the presented qualitative reasoning. 

To complete the zero-amplitude reduction of the modulation system we need to pass to the same limit as a — > in the 
"number of waves" conservation law (|48| in which we assume the aforementioned change of variables (ei, e2, e^, c) 1— > 
(p, u, k, a), 

k z + (uj(p,u, k,a)) x = 0, u = kc. (54) 

As a result, we get 

k z + (uq(p, u, k)) x = 0, (55) 

where 



uj (p,u,k)) = k 5+4/ 7: (56) 




is the dispersion relation (|13jl of linear waves propagating about slowly varying background with locally constant 
values of p and u (here we restrict ourselves with right-propagating waves). Equations (f53"|) . (|55)l comprise a closed 
system which represents an exact zero-amplitude reduction of the full Whitham system (f47|) . (|48| (see [22j], [24[ for 
a detailed justification of this reduction for a class of weakly dispersive nonlinear systems) and, as we shall see, its 
analysis with an account of boundary conditions (fBTil) , ([5T|) yields the necessary information about the leading edge 
x = x + (z) of the dispersive shock wave. 



12 



Now we observe that the "ideal" hydrodynamic equations (|53[) are decoupled from (j55[) and thus, can be solved 
independently for ~p(x, z), u(x, z). However, since the values of p and u at a = are subject to boundary conditions 
(|50| , one should take into account the restriction on the admissible values of p and u at the boundaries of dispersive 
shock wave imposed by the simple- wave transition condition p0[) . Since this restriction is consistent with the equations 
(|53p . it can be incorporated directly into the reduced modulation system by putting 



u = ( arctan vTp — arctan J^y 



(57) 



Substitution of ([57)) into system ()55[) further reduces it to only two differential equations 



p z + v+(p)p x = o, k z + (n(p,k)) x = o, 



(58) 



where 



V + (p) = -^-(arctan \pip~ — arctan */t) H — ^' 
-1/7 1 



7P ' 



(59) 



ft) = w (/9, k) = k 



_V7 V 

The system (|58"|) has two families of characteristics: 




— — ('arctan — arctan ^7) + aI — — y + — . (60) 



|-V +W (61) 



and 

dx dVl{jj, k) 
dz dk 



(62) 



The family (|61[) is completely determined by the simple-wave evolution of the function ~p{x, z) according to the 
dispersionless approximation of the GNLS equation. This family transfers "external" hydrodynamic data into the 
dispersive shock wave region and does not depend on the oscillatory structure. Contrastingly, the behaviour of the 
characteristics belonging to the family (|62|) depends on both p and k. Comparison of the definition (|52|1 of the 
leading edge x + (z) with (|6"2"|) with the account of (|6"0"|) shows that the leading edge of the dispersive shock wave 
represents a characteristic belonging to the family (f6"2"| . Now, since the system ([55)) consists of two equations, then 
according to general prop erties of characteristics of nonlinear hyperbolic systems of partial differential equations (see, 
for instance, one cannot specify two values k and p independently on one characteristic, so the admissible 

combinations of ~p and k at the leading edge of the dispersive shock wave are determined by a characteristic integral 
of the reduced modulation system ([55)) . 

To this end, we substitute k = k(~p) into (J58j to obtain at once 

dk dSl/dp dx dfl 

fl = 0: d= P = v + an/dk on Tz = m- (63) 

The above ordinary differential equation for k must be solved with the initial condition k(p~) = 0. Indeed, since 
the equation (|63() was derived for the case a = it must remain valid in the case of the dispersive shock wave of 
zero intensity, so the dependence k(p) should correctly reproduce the zero wavenumber condition at the trailing edge 
where p = p~ (see ([5T)) ). 
By introducing the variable 



a + (64) 



instead of k, in (|63p . and using Eq. (|60j) . we arrive at the ordinary differential equation 

da (1 + a)[l + 37p + 2a(l - 7p)] 
dp ~ 2p(l+7P)(l + 2a) ' 



(65) 
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with the initial condition 



a(p ) = 1 , 



(66) 



where p is determined in terms of the initial discontinuity (I12p by Eq. (|27p . Once the solution a(p) is found, the 
wave number k + at the leading edge, where ~p = p + = 1, is determined from (|64[) as 



!+= * (1)= «n. 

1 + 7 



(67) 



The velocity of propagation of the leading edge is defined by the kinematic condition (|52[> . which , with an account 
of |62|) . assumes the form 



dx+ _ on 

dz dk 



(l,fc+) 



1 + 7 



2a(l) 



a(l] 



(68) 



For the NLS equation case, i.e. when 7 = 0, the expression for s + in terms of the density jump across the dispersive 
shock wave can be obtained explicitly: the equation 



da 1 



dp 2p 



is readily integrated to give 



<*(?) = 2\ ^-1 
p 



(69) 



(70) 



and thus 



2yV - 1 



for 7 = 



(71) 



in agreement with known results [13j. 

For small values of the saturation parameter 7 <C 1 one can find the correction to this formula with the use of 
Eqs. ([63]) and (|6"5|) . Indeed, if we introduce a — ao + a\, where ao is given by Eq. ([70]) and ai has the order of 
magnitude of 7, then the series expansion of Eq. (|65[) yields the differential equation for the correction a%: 



da 1 
dp 



1 . . ry 5 

27j A^p-fp- If P 
which can be easily solved with account of the initial condition ai(p~) = to give 

4-v/p~-l 1 - Vp~ 1 _ P 



ai(l) = 2 7 Vp^ <^ 1 -p" +64 



hi- 



Then substitution of this expression into Eq. (|68p gives an explicit approximate formula for s -1 



- 1 

which is correct for small 7 as long as ai(l) <C 1. 



-(l- 7 ) + 



(VP" - !) 2 



01 (1) 



(72) 



(73) 



(74) 



2. Trailing edge 

In the vicinity of the trailing edge x = x~(z) the photorefractive dispersive shock wave represents a sequence of 
weakly interacting dark solitons propagating on the slowly varying background p, u. Since one has k — > as x — > x~ , 
we shall be interested in passing to a soliton limit in the modulation system ([4T]) , (|4"B"|) . Instead of performing this 
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limiting passage by a direct calculation (which can be quite involved technically), we shall invoke the reasoning similar 
to that used in the study of the zero-amplitude regime to investigate a reduced modulation system as k — ► 0. 

In limit as k — > 0, the distance between solitons (i.e., a wavelength 2w/k) tends to infinity, so the contribution 
of solitons into the averaged flow p, u vanishes, and, similarly to the case of the vanishing amplitude, we have 
F(p, u) = F(p,u). Hence, we arrive again at the ideal hydrodynamics system (|53[) for p, u. Next, using the arguments 
identical to those used earlier for the case a = but applied now to the case k = we conclude that, for the 
matching condition (|51|) at the trailing edge to be consistent with the simple- wave transition condition (|30p we should 
incorporate the relation (I57p into the reduced as k — > modulation system to obtain the same equation for p (see 
(|58[)). which we reproduce one more time: 

p z + V+(p)p x =0. (75) 

Now we need to pass to the limit as k — > the wave conservation law. This limiting transition, unlike that as 
a — > 0, is a singular one, so it requires a more careful consideration. First we note that the wave conservation law is 
identically satisfied for k = so we need to take into account higher order terms in the expansion of (|54|) for small k. 
Following [22], [HI we introduce a "conjugate wave number" 

/■■=;.-[ / - dp ) (76) 

instead of the amplitude a and the ratio A = k/k instead of the original wave number k, so that the parameters 
(p, u, A, k) provide a new set of the modulation parameters which is convenient for consideration of the vicinity of 
the soliton edge of a dispersive shock. The variable k can be considered as a wave number of oscillations of the variable 
p in the interval e2 < p < £■$ governed by the "conjugate" traveling wave equation 



dp x 2 



= -Q(P) , (77) 



where Q(p) is defined in Eq. (|34[) and 6 is a new phase variable. In the soliton limit e2 — > we can expand Q(p) in 
the vicinity of its minimum point p = ei — &3 so that Eq. (I77|) takes the form of the "energy conservation law" of the 
harmonic oscillator, 



1 d P y 1 d 2 Q 



2 \d6j 4 v 
Then comparison with Eq. (|4ip shows that in this limit 



_(p-p) 2 = Q(p)- 
p 




(78) 

which explains the physical meaning of the variable k in the limit we are interested in. This analogy can be amplified 
by noticing that Eq. ([77]) can be viewed as the traveling wave equation corresponding to the "conjugate" GNLS 
equation obtained from ^ by replacing the variables x and z by ix and iz respectively so that 9 in (|34p is replaced 
by %9 what leads to the change of sign in ([34]) transforming this equation to Eq. ((77)) . Now, the same transformation 
maps a harmonic wave exp[i(fca; — uz)] to the tails of the soliton solution exp[±>t(a; — c so iz)], that is in the soliton 
limit the conjugate frequency ujq can be obtained from the harmonic dispersion relation by a substitution 

iu>o = u>q(ik). (79) 

Actually, this fact is well known and can be used for calculation of the dependence of the soliton velocity c so i = loq/k 
on its inverse half- width k from the dispersion relation for linear waves (see, e.g., [3ll|). Thus, for photorefractive dark 
solitons propagating along the slowly varying background p, u we have the conjugate dispersion relation 
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which, after substitution of the simple-wave relation (f57|) . assumes the form (cf. ([60]) ) 



n(p, k) = u (p, u(p),k) = K 



( arctan 

x/7 v 



7P — arctan 




(81) 



Now we are ready to study the asymptotic expansion as k 
k = Ak into Eq. (J3U) to obtain 



of the wave conservation law (f5"4"| . First we substitute 



kk z + uA x + A{k z + Q x ) = . 



(82) 



where lj = ck. Next we consider Eq. (|8"2")) for small A<1 and assume that A <C A z , A^ for the solutions of our interest 
(this is known to be the case modulation solutions describing dispersive shock waves in weakly dispersive systems, 
where at the soliton edge one has k — > but \k x \, \k z \ — > oo — see [22[ for the general discussion of this behaviour 
and [27j for the detailed calculations in the KdV case). Then to leading order we get the characteristic equation 



OA 

dz 



ndA 

k dx 



= 0. 



which is to say 



A = A 



dx 
dz 



Sl(p, k) 



(83) 



(84) 



where Ao -C 1 is a constant. In particular, when Ao = the characteristic (|84|) specifies the trailing edge (see (|52|) ). 
Now, considering Eq. (|82p along the characteristic family dx/dz = fi/re and using k — n, ui = fl to leading order, we 
obtain 



Kg + fi a 



dx 
dz 



Sl(p, k) 



(85) 



We note that equation n z + Q x = arises as a "soliton wave number" conservation law in the traditional perturbation 
theory for a single soliton (see, for instance, [32]) but to be consistent with the full modulation theory it should be 
considered along the soliton path dx/dz = c so ; =Q,/k. 

Since p and n cannot be specified independently on one characteristic, there should exist a local relationship k(j}) 
consistent with the system (|75p . (1551 . Substituting k — k(~p) into (f85|) and using (|75|) we obtain 



dh~ 



dn/dp 



dp V+- dVL/dn 



(86) 



The initial condition for the ordinary differential equation (|86j) follows from the requirement that the obtained depen- 
dence k(p) should be applicable to the case of the zero-intensity dispersive shock wave, which corresponds to initial 
values p~ = p + = 1. In this case, the width of solitons gets infinitely large, that is k — ► in the limit p — > p + ; this 
follows also from Eq. (|41[) in the limit p m — * pi,. Hence we require k(1) = 0. 

According to the kinematic condition (|52[) the velocity of the soliton edge is equal to the soliton velocity, so we have 



dx 
dz 



Q(p ,K ) 



(87) 



where k = n(p ). 

By introducing a new variable 



s= 4/1- 



k 2 (1 + 7p) 2 



475 



instead of k, Eq. 



reduces to the ordinary differential equation 

da _ (1 + 2) [1 + 37/5+25(1 -7p)] 
dp ~ 2p(l + 7P)(l + 2a) 



(88) 



(89) 
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with the initial condition 

5(1) = 1. (90) 
When the function a{p) is found, the velocity of the trailing soliton is determined by Eqs. (f87|) . (j8lj) . (188|) as 



(arctan xZ-yp — arctan H v ^ _ a(p 

1 + 7P 



(91) 



Then the amplitude a = p — p m of the trailing soliton as a function of the intensity jump p across the dispersive 
shock can be found from the equation (|39p with c = s~, w& = u~ , p/, = 



p-5 2 (p-) _ 2(p--a) 



ja 



1+jp- 



1 



— In 

7a 1 + 7(p — a) 1 + 7P 



(92) 



Again, in the case 7 = corresponding to the NLS equation, all the formulae can be written down explicitly: 
Eq. ([591) reduces to 



(93) 



da 1 + a 
dp 2p 

and its solution satisfying the boundary condition (|9"0"1) is 



a(p) = — - 1. 



Then Eqs. ([HID and dH]) give 



and 



(94) 



(95) 



« = 4(vV-l) 



(96) 



respectively, in agreement with known results [l3j |. Again for small 7 we can find the correction to Eq. (|95| in an 
explicit form. If we denote a = So + Si, where So is given by Eq. (|94|) . then Si satisfies the equation 



dai 
dp 



a x , 8/y/f>- 6 7 



2p y^p-iy/p' 



Si(l) = 0, 



(97) 
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which is readily integrated to give 



) = 



27 



p~ - 1 + 64 



(98) 



and hence 



= VP (l + ai(p )) - 



7- 



(99) 



It is worth noticing that this perturbation approach breaks down for p~ > 16 because of logarithmic divergence 
in Eq. ([98]) as p~ — > 16 — 0. The velocities of the dispersive shock edges as functions of the saturation parameter 
7 are shown in Fig 6. As we see, the presence of even small values of the saturation parameters 7 change the 
expansion velocities considerably compared with the NLS case 7 = because the saturation effects diminish the 
effective nonlinearity which forces the intensive light beam to expand. 



3. Characteristic velocity ordering 



From general point of view, it is important to note that a simple- wave dispersive shock considered above is subject 
to the conditions similar to "entropy" conditions in viscous shocks theory [22, H3 • Basically, these conditions require 
that the number of independent parameters characterizing the modulation solution for the dispersive shock must be 
equal to the number of characteristics families transferring initial data from the x axis into the dispersive shock region 
in the [x, t) plane. For the photorefractive dispersive shock we have four parameters characterizing the initial step 
(fT2")l and one algebraic restriction due to the simple- wave transition condition Thus, the number of independent 
parameters is three. Then, analysis of the characteristic directions at the edges of the dispersive shock waves leads 
to the following inequalities establishing the ordering between the velocities of the dispersive shock edges and the 
characteristic velocities (fl7|) of the dispersionless system : 

VI<8~<V+, V+<s + , s + >s~, (100) 

where subscripts correspond to definitions l|17p and superscripts to two edges of the dispersive shock with constant 
values of and u ± . Inequalities (|100[) provide consistency of the above analytical construction for the derivation of 
the dispersive shock edges, which heavily relies on the properties of characteristics. We have checked that inequalities 
()100p are satisfied for a wide range of parameters. As an illustration, we present in Fig. 7 the plots of the characteristic 
speeds in the simple-wave photorefractive dispersive shock for 7 = 0.2 as functions of the intensity jump across the 
shock. One can see that the ordering (|100p is satisfied. 



18 




FIG. 8: Dependence of the critical intensity p at the trailing edge of the dispersive shock on the saturation parameter 7 for 
fixed value of the intensity p + = 1 at the leading edge. 



a 




FIG. 9: Dependence of the variable 5 on the intensity p~ at the trailing edge of the dispersive shock for fixed value of the 
intensity p + = 1 at the leading edge at 7 = 0.2; the "termination" point corresponds to the intensity p^ r = 4.873 where 
analytical theory loses its applicability. 

4- Vacuum point 

We now investigate dependence of the main properties of the dispersive shock wave on the value of the intensity 
jump across the shock, which is equal to the value p~ at the trailing edge as the value p + = 1 at the leading edge is 
fixed (of course, we assume u + = and u~ given by Eq. (|28|l ). 

It is clear already from the simplest case 7 = that there is a possibility for the value p m at the minimum of the 
trailing dark soliton to become zero (or, which is the same, a = p~) for a certain value of the initial jump p~ . Then 
it follows from (|96[) that this happens at p~ = 4. This gives rise to a vacuum point with p — at the trailing edge 
of the dispersive shock [14| . When the initial step p~ > 4, the vacuum point occurs at some x v inside the dispersive 
shock zone, x~ < x v < x + , and the typical profile of the shock changes (see [HI])- The appearance of the vacuum 
point in the dispersive shock is manifested by the singularity in the profile of u at x = x v but the "momentum" pu 
remains finite. 

For the photorefractive case, when 7 7^ 0, the critical value of p~ = p~ r corresponding to the appearance of the 
vacuum point at the trailing edge of the dispersive shock can be found by putting p~ = a in which immediately 
yields the equation for p~ r 

a(p-) = 0, (101) 

where a(p) is the solution of the ordinary differential equation ([55)1 . The dependence Peril) is shown in Fig. 8. 
Comparison of Eq. (p?Tj) with Eq. <P5|) shows that at the critical point p~ = p~ r we have s _ = u~, that is the trailing 
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soliton is at rest in the reference frame of the intermediate constant state in the decay of an initial discontinuity fT2|) 
(see Section III. A). 

The dependence of a on p~ is shown in Fig. 9. One should note that the change of sign of a at p~ = p cr does 
not constitute nonphysical behaviour even though a as defined by (|88[) is a positive value. In fact, for p~ > p cri 
the velocity u changes its sign at X — Xqj SO that the trailing edge of such a "supercritical" dispersive shock wave 
propagates to the left relative to the vacuum point. To incorporate this change, one should use another branch in the 
linear dispersion relation (| 1 3[) which leads to the change of the sign in the definition of a. As a result, the consistent 
change of signs in (fBTj) and ([59")) leads to the same result for the trailing edge speed s~ defined by (f9"9l . 

One can also see from (1891) that a singularity in the behaviour of a(p~~) is expected at some "termination point" 
p~ = p^ r satisfying 2a(p~) + 1 = for 7 ^ 0. For 7 = 0.2 the value p^ r « 4.873. This singularity has also its 
counterpart in the perturbation theory represented by Eq. (19"5|) . The described pathology in the modulation solution 
for p~ > p^, however, is not confirmed by the direct numerical solutions (see Section IV below) and does not seem 
to have physical sense. One of the explanations of such a discrepancy is that for large values of p~ the accepted 
assumption of applicability of the single-phase modulation theory can fail. Indeed, the developed theory is based 
on the supposition that solutions of our non-integrable photo-refractive system (|11|) behave qualitatively similar to 
their counterparts in the integrable NLS equation case so that the dispersive shock wave can be described with high 
accuracy by the single-phase modulated solution. However, such a supposition can fail in the regions where a drastic 
change of the behavior of a modulated wave takes place. Just this situation occurs in the vicinity of the vacuum 
point, at which the profile of u(x) has a singularity. So one can expect some discrepancy between predictions of the 
modulation theory and exact numerical solutions for the dispersive shock waves with p~ sufficiently close to or grater 
than p cr . As a rough estimate for p~ r one can use the value p~ r = 4 obtained for the integrable NLS equation. Since, 
by definition, ct{Pb r ) = — 1/2 < for all 7 > 0, on can conclude that one always has p7 > p~ r so the predictions of 
the developed modulation theory can become unreliable for such large intensity jumps across the dispersive shock. 



D. Number of solitons generated from a localized initial pulse 

Now we consider an asymptotic evolution of a large-scale decaying initial disturbance 

p(x, 0) = po(x) < 1, u(x, 0) — u a (x); p (x) ^ 1 , u Q (x) — > as |x| — > 00, (102) 

so that the typical spatial scale of this disturbance L> 1. As the numerical simulations for the GNLS equation show, 
such an initial "well" generally decays as z — > 00 into two groups of dark solitons propagating in opposite directions, 
which is consistent with the "two-wave" nature of the GNLS equation. For 7 = the dynamics is described by the 
integrable NLS equation and the soliton parameters are found from the generalized Bohr-Sommcrfeld rule [18j |. In 
the present non-integrable case of the GNLS equation (TlT|) these parameters can be obtained by an extension of the 
modulation method of obtaining the parameters of the dispersive shock wave for the case when the initial distribution 
corresponds to the simple wave solution of the dispersionless equations, that is one of the Riemann invariants (|15|) 
is supposed to be constant. This extension has been developed in [33| in the context of fully nonlinear shallow- 
water waves and we shall use it here to derive the formula for the total number of solitons resulting from the initial 
disturbance l|102[) . First, we assume that for the large-scale initial data (|102p one can neglect the contribution of 
the radiation into the asymptotic as z — > 00 solution, which implies that the whole initial disturbance eventually 
transforms into solitons (this is known to be the case for the integrable NLS equation and is also confirmed by our 
numerical simulations for the GNLS equation). Next, we notice that this transformation into solitons occurs via 
an intermediate stage of the dispersive shock wave formation so we can apply the general modulation theory to its 
description and then to make some inferences pertaining to the eventual soliton train state as z — > 00. 

For definiteness, we consider here the right-propagating dispersive shock wave forming from the profile (|102p satis- 
fying an additional simple-wave restriction 



2 

uo(x) = —— (arctan y/jpo (x) — arctan V7) . (103) 

We consider the wave number conservation law (|54p , which is one of the modulation equations describing the dispersive 
shock wave. For the considered case with decaying at infinity initial profile (|102p we have k — > as |x| — > 00 and, 
therefore, equation l|54|) implies conservation of the total number of waves, 



1 



N = — J kdx = constant. (104) 
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We use an approximate equality sign here due to asymptotic character of the modulation theory which inherently 
cannot predict an integer N exactly. In the Whitham description of the dispersive shock wave, the x-axis is subdivided, 
after the wave breaking at z > Zb, into three regions described in Section III.C : 

— oo < x < x~ (z), x~ (z) < x < x + (z), x + (z) < x < oo , (105) 

where x^(z) are the boundaries of the dispersive shock wave. Generally, these boundaries are not straight lines as in 
the case of the decay of the initial step- like pulse considered above but their nature as characteristics of the modulation 
system remains unchanged. In view of (|105p . the integral in (|104p can be expressed as a sum of three integrals, 

x~ (z) /•x ' {z) />oo I 

k(x, z)dx + / k(x, z)dx + / k(x, z)dx > . (106) 

-oo J x~ (z) J x+ (z) I 

To apply formula (|106p we need first to define the wavenumber k outside the dispersive shock wave as it has been 
actually defined so far only within the nonlinear modulated wave region [x~ (z), x + {z)\. The extended definition of k 
should be consistent with the matching conditions (f50|) . (f5Tj) for all z. 

We know that at the soliton edge x~(z) of the Whitham zone we have k(x~ (z), z) = 0, so we can safely put 
k(x, z) = in the region x < x~(z) and, hence, the first integral vanishes. At the same time, the value of k is 
not explicitly prescribed at the leading edge x + (z) by the boundary condition (|5D|) but is rather determined as a 
function of p due to the fact that the leading edge is a characteristic of the modulation system — see Section III.C. 
The dependence k + (p) is determined then by the ordinary differential equation (|63p (we note that the simple- wave 
transition condition l|28p is already embedded in (f6"3")l and is consistent with the initial conditions (|103p ). This ordinary 
differential equation should be, again, solved with the initial condition k(p~) = 0, and now p~ = pb = 1 where we 
have taken into account that for large pulse (|102[) the wave breaking occurs close to the background intensity, pb = 1. 
Thus, we get the characteristic integral k = k + {p) along the leading edge. The intensity p{x,z) in the downstream 
region x > x + (z) satisfies the simple- wave dispersionless equation 

p z + V+{p)p x = (107) 

with the initial condition p(x,0) = po(x), i.e. the solution is given implicitly by p = p<j(x — V+(p)z). Therefore, to be 
consistent with the boundary values of k prescribed by the characteristic integral of the modulation equations, we have 
to define the wave number downstream the dispersive shock wave as k + (p(x, z)), where p(x,z) is the aforementioned 
simple-wave solution. Then, at z = we get an effective initial distribution of k in terms of the initial data for p 
given by ([ID2J) : 

k(x,0) = k+(po(x)) (108) 

for x > Xb, where Xb is the coordinate of the breaking point; obviously Xb — x~(zb) = x + (zb). Note that this 
definition is also consistent with our definition k = upstream of the dispersive shock wave, since k + (l) = 0. Thus, 
(llOSp describes initial data for the wave number for all x. The function k(x, 0) can be interpreted as the wave number 
distribution for a "virtual" linear modulated wave which accompanies the initial hydrodynamic distributions (p(x,0), 
u(x,0) and transforms, after the wave breaking, into the dispersive shock and, eventually, into a train of solitons. 

Now, we consider (|106p for z — and notice that, since the second integral disappears for z < Zb, Zb > (there is 
no dispersive shock before the breaking point formation so we put x + (z) — x~ (z) for z < Zb), this expression reduces 
to 

k(x,0)dx= / k+(p Q (x))dx. (109) 

-oo J —oo 

As was shown in Section III.C, it is convenient to introduce an auxiliary function a(p) instead of k + (p) according 
to IJMJ) so that 



k + = 2 ^LL± L k (HO) 

1 + IP 

Then a(p) satisfies the ordinary differential equation (|65p with the initial condition a(l) = 1. As a result, the number 
of solitons as z — > oo is determined by the formula 
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FIG. 10: Evolution of the initial step-like pulse with po = 5 and p — 1 for the case of 7 = 0.1. The general structure confirms 
formation of a rarefaction wave, a dispersive shock and an intermediate constant state in between. Intensity p~ = 2.466 
calculated according to Eq. (|27[1 coincides with the numerical result for the intensity of the intermediate state. Coordinates of 
the edges of the rarefaction wave at t = 32 calculated analytically are equal to x^ — —47.7, xf — —9.02 for the rarefaction wave 
and X2 = 42.57, xf = 99.52. One can see that they agree quite well with numerical results. Small-amplitude waves generated 
at around x = — 50 correspond to the linear dispersive "resolution" of the weak discontinuity occurring at the trailing edge of 
the rarefaction wave. 



where ao(x) — a(po(x)). 

When 7 = 0, the solution a(p) of ([65]) is given by ([70]) and assumes here the form 



2 

a = 1 for 7 = 0. (112) 

VP 



Then, for the total number of solitons we have from (jllip 

; o />+oo 



-/ \J{2 - 7^)) 2 - p (x) = - / ^l-pl /2 dx for 7 = 0, (113) 



N 

which agrees with the "simple-wave" reduction of the semi-classical quantization results for the defocusing NLS 
equation obtained in fl8j |. 



IV. NUMERICAL SIMULATIONS OF NONLINEAR WAVES IN PHOTOREFRACTIVE MEDIA 

In this Section, we compare the analytical predictions of the preceding Sections with the results of direct numerical 
simulation of the formation of dispersive shock waves in photorefractive equation ((4]). 

First, we have studied numerically evolution of the step-like pulse. The corresponding results are shown in Fig. 10. 
As we see, all the parameters (velocities of the edges of the rarefaction wave and the dispersive shock, intensity of the 
intermediate state) are in a good agreement with the analytical predictions of Section III. A. 

We have constructed the dependence of p~ and u~ on the saturation parameter 7 using the results of the numerical 
simulations. The results shown in Fig. 11 agree very well with the analytical predictions based on the "simple- wave" 
jump condition (|30[) which is applicable for not too large values of p~ (< 4) such that the vacuum point is not formed. 
In Fig. 12 we show the dependence of the edge "velocities" on the intermediate intensity p~ (with u~ calculated 
according to "simple-wave" jump condition ([28)) ). As we see, a good agreement is observed for p~ < 4. However, 
as p~ increases with growth of po and becomes greater than p~ ~ 4, Eq. (|28|) no longer yields the values of u~ 
compatible with the prescribed value of p~ so that only a single right-propagating dispersive shock is generated; this 
is illustrated by Fig. 13, where a new "intermediate" region of constant flow is seen to be formed which matches with 
the dispersive shock propagating to the right, while another dispersive shock is apparently forming to the left of this 
new constant state providing matching with po . Surprisingly, we have found that the large-amplitude dispersive shock 
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FIG. 11: Dependence of intermediate values (solid lines) of intensity (a) and transverse wave vector (b) on the saturation 
parameter 7 for fixed values of the initial discontinuity parameters: po = 5, uo = for x < and p + = 1, u + = for x > at 
z — 0. Numerically calculated values are shown by crosses. 




FIG. 12: Dependence of s on p (with u calculated according to (J2E|))i 7 — 0.2. Solid lines correspond to analytical formulae 
(|68|) and (|91[) . and dots correspond to results of numeric simulations. 



wave transition between the new intermediate constant state and p = 1 now satisfies a classical shock jump condition 
which follows from the balance of "mass" and "momentum" across the shock as it takes place in classical dissipative 
shocks. Using the dispersionless equations (|14[) represented in a conservative form we find that formal shock jump 
conditions yield the dependence 

- V (n4) 



v /(p- + l)(l+7P)(l+7) 

We have checked that dependence (I114j) is indeed satisfied very well for p~ > 4. The physical mechanism supporting 
the appearance of the classical shock conditions in a dissipationless system such as (|11|) is not quite clear at the 
moment. We note that similar effect of appearance of the classical shock jump condition across the expanding 
dispersive shock have been recently observed in (34j for large-amplitude shallow-water undular bores modeled by the 
Green-Naghdi system, which is also not integrable by the 1ST. At the same time, it is known very well that for the 
dispersive shocks described by the integrable NLS equation, the simple-wave jump condition is satisfied exactly for all 
values of initial density jump - this follows from the full modulation solution (see [HI, [2(|, Q) and is also confirmed 
by our numerical simulations. So it is possible that the described phenomenon of the appearance of the classical shock 
conditions constitutes a specific manifestation of nonintegrability in dispersive dissipationless systems which is yet to 
be explored. 

Next, we have compared the analytical predictions of subsection III.D for number of dark solitons generated from 
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FIG. 13: Dispersive shock evolving from the step-like pulse with p~ and u~ related by the "simple-wave" jump condition 
for large value of p~ = 10 much greater than p~ = 4. Occurrence of a vacuum point in the region between x = 100 and 
x — 150 is clearly seen. A new intermediate constant state is formed in the region behind the dispersive shock showing that 
the simple- wave jump condition (|28[1 does not prevent anymore the formation of the second wave for large values of p~ . 



t 

t 



■■ 

.100 



r 



-200 



-100 



100 



200 



FIG. 14: Profile of intensity at z = 100 evolved from the initial pulse (|115|l (dashed line) with initial profile of u(x) calculated 
according to (|103[) . 



a hole-like disturbance with numerical simulations. We took the initial distribution of intensity 

^^{'-^k^) 2 (115) 

and the initial distribution of transverse wave number was calculated according to the equation (|103p . Evolution 
of such a pulse according to the photorefractive equation with 7 = 0.2 is illustrated by Fig. 14 where the profile of 
intensity is shown at z = 100. As we see, this pulse, after the wave breaking and formation of a dispersive shock, 
evolves eventually into a number of dark solitons. We note that the appearance of several solitons propagating to 
the left does not contradict to the unidirectional restriction guaranteed by the simple- wave initial conditions (|115[) . 
(|103[) - these left-propagating solitons occur due to relatively high amplitude of the initial disturbance (|115[) . which 
leads to the appearance of the vacuum point at the intermediate stage of the dispersive shock wave and, therefore, to 
the formation of some number of left-propagating solitons - see Section III.C. The total number of created solitons 
calculated by means of the modulation formula (jllip as a function of the saturation parameter 7 is shown by solid 
line in Fig. 15 and the corresponding results of numerical simulations are indicated by dots. Taking into account the 
asymptotic nature of the developed analytical theory for this integer-valued function, and the fact that considered 
initial data (| 1 1 5[) produce a vacuum point (i.e. at some stage of the dispersive shock development the "instantaneous" 
initial jump p~ > p~ r ), the agreement can be considered as quite good. 
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FIG. 15: Number of solitons N as a function of 7; solid line represents analytical dependence and dots correspond to 

numerical simulations. 
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FIG. 16: Density plot of the output intensity evolved from a strip-like initial distribution. 



In Refs. p, Q the theory of dispersive shocks, evolving from a step-like pulse according to the NLS equation ([5]) 
(7 = 0) was used for qualitative explanation of dispersive shocks with other geometries in concrete physical situations 
(see also [Il[ where the NLS theory of the wave breaking was also used for description of dispersive shocks in Bose- 
Einstcin condensates). In a similar way, the developed here theory of dispersive shocks in photorefractive media can 
be used for the description of experiments on generation of optical shocks. Such experiments were described in @, Q 
and here we present some results based on the numerical solutions of the photorefractive equation (j4|) with the initial 
conditions similar to the initial light distributions in the mentioned experimental works (similar results of numerical 
simulation were presented in Q ) . 

In [8| the distribution of output intensities are presented for initial distributions in the form of a strip, a circle, 
and two separated circles. We have performed numerical simulations with similar initial conditions. In Fig. 16 we 
present a density plot of the output intensity evolved, according to the photorefractive equation with 7 = 0.1, from 
the strip-like initial distribution given by the formula 



p(x) = 



1 + 5(1 
1 



x 2 /25) - 2 for 
for 



\x\ < 5, 
\x\ > 5, 



(116) 



which approximates with a good enough accuracy the constant values of intensities inside the strip and outside it. 
Similar density plot for the circle initial distribution is shown in Fig. 17. 

As we see, in both cases the initial "hump" breaks with formation of dispersive shocks — in the strip-like geometry 
we get two shocks propagating in opposite directions and in circular geometry we have a ring-like dispersive shock 
expanding in radial direction. 
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FIG. 17: Density plot of the output intensity evolved from a circle initial distribution. 
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FIG. 18: Interaction of two circular dispersive shocks. 



In Fig. 18 an interaction of two circular dispersive shock waves is shown. It is remarkable that even in this two- 
dimensional nonintegrable photorefractive case, the nonlinear dispersive shock waves interact apparently elastically, 
without production of other waves in the region of their overlap. It is this kind of picture that is expected in the 
system with 7 = described by the integrable NLS equation where the interaction of two dispersive shocks leads to the 
formation of a two-phase modulated wave region described by the corresponding multiphased-averaged modulation 
system [l{|. While the analytical description of multiphase nonlinear waves in photorefractive equation (jTTJ) is not 
available, the qualitative similarity between the solution behaviour for the nonintegrable photorefractive equation and 
for the NLS equation for moderate values of initial amplitudes can be considered as a confirmation of robustness of 
the modulated travelling wave ansatz in the description of dispersive shock waves in nonintegrable systems, at least 
for some reasonable range of initial amplitudes. 



V. CONCLUSION 



In this paper, we have developed the theory of formation of dispersive shocks in propagation of intensive light beams 
in photorefractive optical systems. The theory is based on Whitham's modulation approach in which a dispersive 
shock is described as a modulated nonlinear periodic wave and slow evolution along the propagation axis is governed 
by the averaged modulation equation. In spite of the absence of complete integrability of the equation describing 
propagation of light beams in photorefractive media, the main characteristic parameters of shocks can be determined 
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by means of the approach developed in [22|-[24( and based on the study of reductions of Whitham equations for the 
wave regimes realized at the boundaries of the dispersive shock. In particular, "velocities" of the dispersive shock 
edges are found as functions of the jump of intensity across the shock as well as amplitude of the soliton at the rear 
edge of the shock. The number of solitons produced from a finite initial disturbance is also determined analytically 
for initial distributions related by a so-called simple wave condition. The analytical theory agrees very well with 
numerical simulations as long as there is no vacuum point in the shock. Appearance of a vacuum point leads to the 
formation of a singularity in a "transverse" wavevector distribution and such a drastic change in the wave behavior 
cannot be traced by the developed approach. However, this situation occurs at very high input intensities of a light 
beam so that for practical purposes the developed theory provides accurate enough approximation. 

Although the theory is essentially one-dimensional (i.e. with one transverse space coordinate) it can give qualitative 
explanation of experiments with other geometries, which is illustrated by the results of numerical simulations. 
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